namespace NBody.InitialConditions
{
    using System;
     
    public class Sphere
    {
        double r_scale;
        double v_scale;
        double total_mass;
        
        bool hubble;
        
        public Sphere(double radius, double v_scale, double total_mass)
        {
            this.r_scale = radius;
            this.v_scale = v_scale;
            this.total_mass = total_mass;
            this.hubble = false;
        }
        
        public double Radius
        {
            get { return r_scale; }
        }
        
        public double Vscale
        {
            get { return v_scale; }
        }
        
        public double TotalMass
        {
            get { return total_mass; }
        }
        
        public bool InHubbleFlow
        {
            get { return hubble; }
            set { hubble = value; }
        }
        
        public NBodySystem GenerateGrid(int n_per_side, double time)
        {
            Cube c = new Cube(2.0 * r_scale, 6.0 / Math.PI * total_mass, Cube.Origin.Centre);
            NBodySystem s = c.GenerateGrid(n_per_side, time);
            
            s.ExtractSphere(r_scale);
            
            if (hubble)
            {
                for (int i = 0; i < s.NumParts; i++)
                {
                    double vx = v_scale * s[i].X;
                    double vy = v_scale * s[i].Y;
                    double vz = v_scale * s[i].Z;
                    
                    s[i].SetVelocity(vx, vy, vz);
                }
            }
            
            return s;
        }
        
        // generate system from distribution function
        // see Aarseth et al for implementation details
        
        public delegate double DFFunction(double x);
        
        public NBodySystem GenerateFromDF(int n, double time, DFFunction mass_function, DFFunction potential_function, DFFunction g_function, int seed)
        {
            NBodySystem s = new NBodySystem(n, time);
            Random random;
            if (seed > 0)
                random = new Random(seed);
            else
                random = new Random();
            
            double m_frac = 0.999; // assumes total_mass = 1 ?
            double g_max = 0.1; // no, this depends on g_function
                                    
            for (int i = 0; i < n; i++)
            {
                double mass = total_mass / (double)n;
                double radius = r_scale * mass_function(m_frac * random.NextDouble());
                
                double x, y, z;
                PickShell(random, radius, out x, out y, out z);
                
                double q, test;
                do
                {
                    q = random.NextDouble();
                    test = g_max * random.NextDouble();
                } while (test > g_function(q));
                
                double velocity = v_scale * q * Math.Sqrt(-2.0 * potential_function(radius));
                
                double vx, vy, vz;
                PickShell(random, velocity, out vx, out vy, out vz);
                
                s.AddParticle(new Particle((double)mass, (double)x, (double)y, (double)z, (double)vx, (double)vy, (double)vz));
            }
            
            return s;
        }
        
        private double plummer_mass(double x)
        {
            return 1.0 / Math.Sqrt(Math.Pow(x, -2.0/3.0) - 1.0);
        }
        
        private double plummer_potential(double x)
        {
            return -1.0 / Math.Sqrt(1.0 + x * x);
        }
        
        private double plummer_g(double x)
        {
            return x * x * Math.Pow(1.0 - x * x, 3.5);
        }        
        
        public NBodySystem GeneratePlummer(int n, double time, int seed)
        {
            return GenerateFromDF(n, time, new DFFunction(plummer_mass), new DFFunction(plummer_potential), new DFFunction(plummer_g), seed);
        }
        
        private void PickShell(Random random, double radius, out double x, out double y, out double z)
        {
            double theta = Math.Acos(1.0 - 2.0 * random.NextDouble());
            double phi = 2.0 * Math.PI * random.NextDouble();
    
            x = radius * Math.Cos(phi) * Math.Sin(theta);
            y = radius * Math.Sin(phi) * Math.Sin(theta);
            z = radius * Math.Cos(theta);
        }
    }
}
                
                
            
            
            
        
            
